#This is to obtain:
#table_A17.tex
#figure_A14a.jpg
#figure_A14b.jpg
#figure_A14c.jpg
#figure_A14d.jpg

rm(list = ls())
library("ggplot2")
library("lfe")
library("readr")
library("sf")
library("lemon")
library("stringi")
library("readxl")
library("dplyr")
library("stargazer")


setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")

#Step1: Bosnia
#Reading Bosnia shape
BIH <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
               layer="BIH_adm0")

BIH_roads <- st_read(dsn="./Dropbox/Legacies_Project/Data/Roads/roads_boasnia.gdb",
                       layer="BIH_adm2")
BIH_roads<-st_simplify(BIH_roads, preserveTopology = FALSE, dTolerance = 0.005)

#Adding and editing variables
BIH_roads$NAME_3=stri_trans_general(BIH_roads$NAME_3, "Latin-ASCII")
BIH_roads$NAME_3<-toupper(BIH_roads$NAME_3)
BIH_roads$treat<-NA
BIH_roads$treat<-0

#Labelling border districts
BIH_roads$treat[BIH_roads$NAME_3=="BIHAC"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BIJELJINA"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BOSANSKA GRADISKA"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BOSANSKA KOSTAJNICA"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BOSANSKA KRUPA"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BOSANSKA KRUPA"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BOSANSKI BROD"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BOSANSKO GRAHOVO"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BRCKO"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="BUZIM"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="CAZIN"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="DERVENTA"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="DOMALJEVAC-SAMAC"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="DRVAR"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="DUBICA"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="NOVI GRAD"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="ODZAK"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="ORASJE"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="SRBAC"]<-1
BIH_roads$treat[BIH_roads$NAME_3=="VELIKA KLADUSA"]<-1

BIH_roads$croatia_border_area<-BIH_roads$treat
BIH_roads$road_density_2022<-BIH_roads$road_density


#Step2: Croatia
x<-st_layers(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb")
HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")
krajna <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)

final_sp_incl_brd<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
final_sp_incl_brd = st_transform(final_sp_incl_brd, st_crs(BIH_roads))
final_sp_incl_brd<-st_simplify(final_sp_incl_brd, preserveTopology = FALSE, dTolerance = 0.005)

#Step3: Analyses
#Sample1
x1<-lm(road_density_2022~treat, data = final_sp_incl_brd)
summary(x1)
mean_cr_mil<-round(mean(final_sp_incl_brd$road_density_2022, na.rm = T),2)
sd_cr_mil<-round(sd(final_sp_incl_brd$road_density_2022, na.rm = T),2)
final_sp_incl_brd$bosnia_border_mun[is.na(final_sp_incl_brd$bosnia_border_mun)]<-1
final_sp_incl_brd$treat[is.na(final_sp_incl_brd$treat)]<-1
x<-subset(final_sp_incl_brd, select = c("treat", "bosnia_border_mun"))


#Sample2
new2<-subset(final_sp_incl_brd, treat==0 | (bosnia_border_mun==1 & treat==1))
new2$treat<-new2$bosnia_border_mun
x2<-lm(road_density_2022~treat, data = new2)
summary(x2)
mean_cr_strip<-round(mean(new2$road_density_2022, na.rm = T),2)
sd_cr_strip<-round(sd(new2$road_density_2022, na.rm = T),2)

#Sample3
BIH_roads$croatia_border_area

x3<-lm(road_density_2022~croatia_border_area, data = BIH_roads)
summary(x3)
mean_bih_strip<-round(mean(BIH_roads$road_density, na.rm = T),2)
sd_bih_strip<-round(sd(BIH_roads$road_density, na.rm = T),2)

#Sample4
BIH_roads2<-subset(BIH_roads, croatia_border_area ==1)
final_sp_incl_brd2<-subset(final_sp_incl_brd, select = c(treat , road_density_2022))
final_sp_incl_brd2<-subset(final_sp_incl_brd2, treat==1)
BIH_roads2<-subset(BIH_roads2, select = c(treat, road_density_2022))
BIH_roads2$treat[BIH_roads2$treat==1]<-0
HRV_BIH<-rbind(x = BIH_roads2, y = final_sp_incl_brd2)

x4<-lm(road_density_2022~treat, data = HRV_BIH)
summary(x4)
mean_hrvbih_strip<-round(mean(HRV_BIH$road_density, na.rm = T),2)
sd_hrvbih_strip<-round(sd(HRV_BIH$road_density, na.rm = T),2)



column_labels= c("\\shortstack{Croatia}", 
                 "\\shortstack{Croatia}",
                 "\\shortstack{Bosnia}",
                 "\\shortstack{Bosnia and\\\\ Croatia}")


roads_placebo =stargazer(x1,
                        x2,
                        x3,
                        x4,
                        column.labels= column_labels,
                        keep = c("treat",
                                 "croatia_border_area"),
                        covariate.labels=c("Military Territory", 
                                           "Border Territory in Bosnia"),
                        dep.var.caption = "Dependent variable: Road Density, 2022
                                Specification:",
                        dep.var.labels.include = F,
                        omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                        add.lines=list(c("Mean", mean_cr_mil, 
                                         mean_cr_strip, 
                                         mean_bih_strip,
                                         mean_hrvbih_strip),
                                       c("SD", sd_cr_mil, 
                                         sd_cr_strip, 
                                         sd_bih_strip,
                                         sd_hrvbih_strip)))
note_latex <- "\\multicolumn{5}{l} {\\parbox[t]{15cm}{ \\footnotesize \\textit{Notes:}
Coefficients and standard errors in parantheses from OLS regression.
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01. The unit of analysis is municipality (opcina).
Data captures modern road density from Open Street Map.}}"
roads_placebo [grepl("Note", roads_placebo)] <- note_latex
#cat(pct_no_water, file="pct_no_water.tex", sep="\n")
cat(roads_placebo, file="./Dropbox/Legacies_Project/Paper/tables/table_A17.tex", sep="\n")

######
#Maps#
######

#Sample1
HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")
HRV_adm0<-st_simplify(HRV_adm0, dTolerance = 0.005)

mil_col_short_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")
mil_col_short_polygon = st_transform(mil_col_short_polygon, st_crs(HRV_adm0))
mil_col_short_polygon<-st_simplify(mil_col_short_polygon,  dTolerance = 0.005)

all_coun<- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                   layer="all_countries_GCS")
all_coun = st_transform(all_coun, st_crs(HRV_adm0))
all_coun<-st_simplify(all_coun, dTolerance = 0.005)

cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))

final_sp_incl_brd_yes<-subset(final_sp_incl_brd, treat==1)
final_sp_incl_brd_no<-subset(final_sp_incl_brd, treat==0)


map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "black", linewidth = 1)+
  geom_sf(data = all_coun, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", size = 0.06)+
  geom_sf(data = final_sp_incl_brd_yes, aes(fill="Military Territory"),
          size = 0.1)+
  geom_sf(data = final_sp_incl_brd_no, aes(fill="Civilian Territory"),
          size = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, size = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  ggtitle("(1) Croatia")+
  labs(x = "Longitude", y="Latitude")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12),
    plot.margin = unit(c(0, 0, 0, 0), "null"),
    panel.spacing = unit(c(0, 0, 0, 0), "cm"))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/figure_A14a.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)


#Sample 2
cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))


final_sp_incl_brd$bosnia_border_mun
final_sp_incl_brd_yes_brd<-subset(final_sp_incl_brd, bosnia_border_mun==1 & treat==1)
final_sp_incl_brd_no_brd<-subset(final_sp_incl_brd, bosnia_border_mun==0 & treat==0)


map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "black", linewidth = 1)+
  geom_sf(data = all_coun, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", size = 0.06)+
  geom_sf(data = final_sp_incl_brd_yes_brd, aes(fill="Military Territory"),
          size = 0.1)+
  geom_sf(data = final_sp_incl_brd_no_brd, aes(fill="Civilian Territory"),
          size = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, size = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  ggtitle("(2) Croatia")+
  labs(x = "Longitude", y="Latitude")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12),
    plot.margin = unit(c(0, 0, 0, 0), "null"),
    panel.spacing = unit(c(0, 0, 0, 0), "cm"))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/figure_A14b.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)


#Sample3
cols.fill = c("Border Territory in Bosnia" = alpha("yellow", 0.5), 
              "Other Territory in Bosnia" = alpha("chartreuse4", 0.5),
              "Excluded" = alpha("gray", 0.01))


BIH_roads_yes_brd<-subset(BIH_roads, croatia_border_area==1)
BIH_roads_no_brd<-subset(BIH_roads, croatia_border_area==0)


map<-ggplot() + 
  geom_sf(data = BIH, fill = NA, colour = "black", linewidth = 1)+
  geom_sf(data = all_coun, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = BIH_roads_yes_brd, aes(fill="Border Territory in Bosnia"),
          size = 0.1)+
  geom_sf(data = BIH_roads_no_brd, aes(fill="Other Territory in Bosnia"),
          size = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  ggtitle("(3) Bosnia")+
  labs(x = "Longitude", y="Latitude")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12),
    plot.margin = unit(c(0, 0, 0, 0), "null"),
    panel.spacing = unit(c(0, 0, 0, 0), "cm"))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/figure_A14c.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)


#Sample 4
cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Border Territory in Bosnia" = alpha("yellow", 0.5),
              "Excluded" = alpha("gray", 0.01))


HRV_BIH_yes<-subset(HRV_BIH, treat==1)
HRV_BIH_no<-subset(HRV_BIH, treat==0)


map<-ggplot() + 
  geom_sf(data = BIH, fill = NA, colour = "black", linewidth = 1)+
  geom_sf(data = HRV_adm0, fill = NA, colour = "black", linewidth = 1)+
  geom_sf(data = all_coun, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = BIH_roads, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = HRV_BIH_yes, aes(fill="Military Territory"),
          size = 0.1)+
  geom_sf(data = HRV_BIH_no, aes(fill="Border Territory in Bosnia"),
          size = 0.1)+
  geom_sf(data = krajna, colour = "blue", fill = NA, size = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  ggtitle("(4) Bosnia and Croatia")+
  labs(x = "Longitude", y="Latitude")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12),
    plot.margin = unit(c(0, 0, 0, 0), "null"),
    panel.spacing = unit(c(0, 0, 0, 0), "cm"))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/figure_A14d.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)
